############################
#
# Create plots
#
############################
library(plotrix)
library(rdrobust)
library(readstata13)

#############################
#
# US Senate
#
#############################

data = read.dta13("data-senate.dta")
names(data)

# RD plot
pdf("./final-rdplot-senate.pdf")
rdplot(data$demwinnext2, data$demmv, x.label="Margin of Victory at t", y.label="Victory at t+1", title="Democratic Party",  x.lim=c(-50,50), y.lim=c(0,1), binselect = "esmv", scale=1, lowerend=-100,upperend=100)
dev.off()

# Histogram
quantile(data$votesh_so[abs(data$demmv)<3], prob=c(0.01,0.05,0.10,0.15,0.20,0.50,0.75,0.90))
p1 = 48.5
q1=mean(data$votesh_so[abs(data$demmv)<3] < p1) * 100
p2 = 46
q2=mean(data$votesh_so[abs(data$demmv)<3] < p2) * 100

pdf("./final-histogram-stropp-senate.pdf")
hist(data$votesh_so[abs(data$demmv)<3], xlab="Vote Share (%)", col="cyan", main="Democratic Party's Strongest Opponent", axes=F, xlim=c(20,52), breaks=7, ylab="")
axis(side=1,at=c(20,25,30,35,40,p1,p2,50,52), labels=paste(c(20,25,30,35,40,p1,p2,50,52)))
axis(side=2,at=seq(0,70,10), labels=paste(seq(0,70,10)))
abline(v=p1, col="blue", lty=2,lwd=2)
arrows(40,50,p1,50,length=0.10)
text(x=36.7,y=50,labels=paste(round(q1,0),"th percentile", sep=""))
abline(v=p2, col="red", lty=2,lwd=2)
arrows(40,30,p2,30,length=0.10)
text(x=36.7,y=30,labels=paste(round(q2,0),"th percentile", sep=""))
dev.off()

#############################
#
# Brazil PSDB 
#
#############################
data = read.dta13("data-election-example.dta")
names(data)
# RD plot
pdf("./final-rdplot-psdb.pdf")
rdplot(data$y, data$x_norm, x.label="Margin of Victory at t", y.label="Victory at t+1", title="PSDB", x.lim=c(-50,50), y.lim=c(0,1), binselect = "esmv", scale=1, p=4)
dev.off()

# Histogram
quantile(data$votesh_so[abs(data$x_norm)<3], prob=c(0.01,0.05,0.10,0.15,0.20,0.50,0.75,0.90))
p1 = 48.5
q1=mean(data$votesh_so[!is.na(data$votesh_so) & abs(data$x_norm)<3] < p1) * 100
p2 = 46
q2=mean(data$votesh_so[!is.na(data$votesh_so) & abs(data$x_norm)<3] < p2) * 100

pdf("./final-histogram-stropp-psdb.pdf")
hist(data$votesh_so[abs(data$x_norm)<3], xlab="Vote Share (%)", col="cyan", main="PSDB's Strongest Opponent", breaks=20, axes=F, xlim=c(20,52))
axis(side=2,at=seq(0,500,100), labels=c(paste(seq(0,500,100))))
axis(side=1,at=c(20,25,30,35,40,p2,p1,52), labels=paste(c(20,25,30,35,40,p2,p1,52)))
abline(v=p1, col="blue", lty=2,lwd=2)
arrows(40,250,p1,250,length=0.10)
text(x=36.8,y=250,labels=paste(round(q1,0),"th percentile", sep=""))
abline(v=p2, col="red", lty=2,lwd=2)
arrows(40,150,p2,150,length=0.10)
text(x=36.8,y=150,labels=paste(round(q2,0),"th percentile", sep=""))
legend(19,400,"48.5 = Minimum value possible in perfect two-party system", col="blue", lty=2,lwd=2, cex=0.9)
dev.off()

###################################################
#
# Create final plots using rdrboust results generated by ./create-postfiles-for-R.do
#
###################################################
data = read.dta13("RD-results.dta")

########################
# Election example
########################
results = data[(data$example == "election" & data$analysis=="cutoff-by-cutoff"),]
y = results[,"tau"]
x = results[,"cutoff"]
taup = results[1,"taup"]
pdf.file = "./final-cutoff-by-cutoff-plot-election.pdf"
xlab = "Vote Share of Strongest Opponent"
ylab =  "RD Effect on Electoral Victory"
gr = results$cutoff
xlim=c(33,50); ylim=c(-1,1); xticks=gr; yticks=seq(-1,1,0.2)

pdf(pdf.file) 
plotCI(x, y, xlab=xlab,ylab=ylab, xlim=xlim, ylim=ylim, axes=FALSE, ui=results[,"CIurb"], li=results[,"CIlrb"])
axis(side=1,at= xticks, labels=paste(xticks, sep=""))
axis(side=2,at= yticks, labels=paste(yticks, sep=""))
abline(h=0, col="gray", lty=1,lwd=2)
abline(h=taup, col="blue", lty=2,lwd=2)
predict(lm(y~x+x^2))
lines(x,predict(lm(y~x+I(x^2))), col="red", lty=2, lwd=1.5)
legend("topleft", legend=c("Pooled estimate", "Quadratic polynomial fitted to mean effects"),
       lty=c(2,2), lwd=c(2,1.5),col=c("blue", "red"))
dev.off()    

########################
# Education example
########################
results = data[(data$example == "education" & data$analysis=="cutoff-by-cutoff"),]
cutoffs = results$cutoff
y = results[,"tau"]
x = results[,"cutoff"]
taup = results[1,"taup"]
yticks = c(-30, -20, -10, 0, 10)
xlab = "Cutoff Value (Region-specific Past Test Scores Index)"
ylab = "RD Effect on Language Test Score"
pdf.file = "./final-cutoff-by-cutoff-plot-education.pdf"
xlim=c(42,52)
ylim=c(-35, 20)
pdf(pdf.file)
plotCI(x, y, xlab=xlab,ylab=ylab, xlim=xlim, ylim=ylim, axes=TRUE, ui=results[,"CIurb"], li=results[,"CIlrb"])
abline(h=0, col="gray", lty=1,lwd=2)
abline(h=taup, col="blue", lty=2,lwd=2)
predict(lm(y~I(x)+I(x^2)))
lines(x,predict(lm(y~I(x)+I(x^2))), col="red", lty=2, lwd=1.5)
legend(x=42.92,y=-20, legend=c("Pooled estimate", "Quadratic polynomial fitted to mean effects"),
       lty=c(2,2), lwd=c(2,1.5),col=c("blue", "red"), cex=0.9)
dev.off()


########################
# Population example
########################
results = data[(data$example == "population" & data$analysis=="cutoff-by-cutoff"),]
y = results[,"tau"]
x = results[,"cutoff"]
taup = results[1,"taup"]
yticks = c(-30, -20, -10, 0, 10)
xlab = "Cutoff Value (Population)"
ylab = "RD Effect on Corruption Indicator"
pdf.file = "./final-cutoff-by-cutoff-plot-population.pdf"
xlim=c(9000,45000)
ylim=c(-2,2)
pdf(pdf.file)
plotCI(x, y, xlab=xlab,ylab=ylab, xlim=xlim, ylim=ylim, axes=TRUE, ui=results[,"CIurb"], li=results[,"CIlrb"])
abline(h=0, col="gray", lty=1,lwd=2)
abline(h=taup, col="blue", lty=2,lwd=2)
predict(lm(y~I(x)+I(x^2)))
lines(x,predict(lm(y~I(x)+I(x^2))), col="red", lty=2, lwd=1.5)
legend("topleft", legend=c("Pooled estimate", "Quadratic polynomial fitted to mean effects"),
       lty=c(2,2), lwd=c(2,1.5),col=c("blue", "red"), cex=0.9)
dev.off()
